Adaptation and Performance of the Cartesian 
Coordinates Fast Multipole Method for Nanomagnetic 

Simulations 

Wen Zhang, Stephan Haas 

Department of Physics and Astronomy, University of Southern California 
Los Angeles, CA 90089, USA 



Abstract 

An implementation of the fast multiple method (FMM) is performed for mag- 
netic systems with long-ranged dipolar interactions. Expansion in spherical 
harmonics of the original FMM is replaced by expansion of polynomials in 
Cartesian coordinates, which is considerably simpler. Under open boundary 
conditions, an expression for multipole moments of point dipoles in a cell is 
derived. These make the program appropriate for nanomagnetic simulations, 
including magnetic nanoparticles and ferrofiuids. The performance is opti- 
mized in terms of cell size and parameter set (expansion order and opening 
angle) and the trade off between computing time and accuracy is quantita- 
tively studied. A rule of thumb is proposed to decide the appropriate average 
number of dipoles in the smallest cells, and an optimal choice of parameter 
set is suggested. Finally, the superiority of Cartesian coordinate FMM is 
demonstrated by comparison to spherical harmonics FMM and FFT. 
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1. Introduction 

Magnetic nanoparticles and fluids have attracted intensive investigation 
during the last two decades [H, 0, 0, 0, 0]. Micromagnetic simulation is an 
powerful tool to study such systems, where the main challenge is the dipolar 
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interaction between magnetic moments, which can be dealt with either dis- 
cretely or in the continuum limit. Since the interaction is long ranged, the 
complexity of a brute-force calculation will be 0(N 2 ). One popular method 
to improve this performance is by Fast Fourier Transform (FFT), which re- 
duces the complexity to 0(N * log(N)). However, FFT suffers from the fact 
that it requires a regular lattice arrangement of dipoles, and a large num- 
ber of padding areas have to be added when dealing with exotic geometries 
and open boundary conditions. Thus, the alternative, fast multipole method 
fFMM] has attracted increasingly more attention in the past several years 

la SB- 

The FMM was first introduced by L. Greengard et al [lOj, LL3J, and has 
been used ever since to speed up large scale simulations involving long ranged 
interactions. It has the charming advantage of O(N) complexity. Further- 
more, there is no constraint on the distribution of particles and the bound- 
aries of the system. Hence it is extremely useful for magnetic fluid systems 
and nanomagnets with exotic geometries. Last but not least, it is a scal- 



able algorithm, i.e. it can be efficiently implemented in parallel [121 ] . With 
so many advantages though, the adoption of FMM in magnetic system is 
not widely applied, probably due to two reasons. One concern is that since 
there is a huge overhead to achieve O(N) complexity, the FMM becomes 
faster only for very large N. Actually this is true only when one considers 
a very high order expansion (say 10). In the context of micromagnetics, we 
will show that an expansion to the order of no more than 6 will be satisfy- 
ing, and with some optimization procedure FMM will be superior even in a 
system of 10 3 dipoles. The other reason is that the standard implementa- 
tion of the FMM algorithm is based on the well-known spherical harmonics 
expansion 13( of 1/r. For dipolar interaction which decays as 1/r 2 , how- 



ever, it is not straightforward to apply. Further, the additional complexity of 
calculating spherical harmonics acts as another barrier. To overcome these 



shortcomings, expansions in Cartesian coordinate were proposed [14|, [15|, but 



they were not applied successfully until recently [8f]. Following these studies, 
here we present a very simple but useful formula (Eq. [3]) to calculate the 
multipole moments in dipolar systems. 

In Section II and III, a brief description of the FMM algorithm and all 
the necessary equations is presented. In Section IV, performance issues will 
be discussed in detail. The number of dipoles contained in the smallest cell 
will be shown to be crucial in terms of performance and a rule of thumb on 
choosing the right number will be given (Eq. [T7|) . Then, the optimal choice 



2 



of expansion order and opening angle to achieve certain error bounds will be 
discussed, followed by quantitative study of the trade off between computing 
time and accuracy. Finally, we compare the Cartesian coordinates FMM 
with spherical harmonics FMM and FFT. 



2. Algorithm 

The crucial ideas of the FMM are: (i) chunk source points together into 
large cells whose field in remote cells is computed by a multipole expansion 
of the source points; (ii) use a single Taylor expansion to express the smooth 
field in a given cell contributed by the multiple expansions of all remote cells. 
As the details of the algorithm are discussed in previous literature 10, 0], 



here we will only briefly describe them, stressing those aspects which are 
special to our implementation. Instead of the regular geometric hierarchi- 
cal traversals of the system, we follow a simple alternative way[?J by using 
recursive functions. This method simplifies the program significantly and is 
especially elegant for system with exotic geometries, since one only needs to 
recursively divide the system into halves. 

After setting up all the cells, the entire program mainly consists of three 
recursive functions traversing the entire system trice. Pseudocodes are pro- 
vided in the Appendix. 

1. Downward Pass: construct "partners" list (similar to the interaction 



list in the original FMM [10J). For each cell, the partners list contain those 
cells which are near to its parent but far to itself. Instead of the original 
geometric rules (near cell pairs are those cells touching each other), near 
cell pairs refer to those with opening angle a larger than certain value a m 
(a m < 1). The opening angle is defined as 2r/d (a m = 2r/d), where r is 
the "radius" of two cells and d is the distance between their centers. Note 
that "radius" denotes the maximum distance between a corner of one cell 
and its center. With these definitions, it is easy to construct the partners 
lists. First, set the partners list of each cell to be empty. Then, for each 
cell except those cells without children (leaf cells), start from its partners list 
generated by its parent (for the root cell, this list is empty) and perform the 
following procedure: (i) add the children of all its near cells it to the partners 
list of each of its children, and delete these cells from its own partners list; 
(ii) add one of its children to the other child's partners list (do it for both 
children); (iii) call its children to perform these tasks. In this way, we start 
from the root cell and recursively traverse its children until reaching leaf 
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cells. For each leaf cell, add to another list called "nearPartners" those cells 
in the partners list which are near to it and delete them from its partners 
list. Finally put itself in its nearPartners list. 

2. Upward Pass: construct the multipole moments for each cell. Start 
from the leaf cells and explicitly calculate their multipole moments according 
to Eq. [3] below. Then upward traverse all the other cells in binary tree 
hierarchy whose multipole moments are obtained by shifting the origin of its 
children's multipole moments (see Eq. [7]). This process can be easily realized 
by a recursive function. 

3. Downward Pass: construct the Taylor expansion of the smooth field in 
each cell. Start from the root cell and down traverse all the other cells in the 
binary tree hierarchy. For each cell, inherit the Taylor expansion coefficients 
of the smooth field from its parents according to Eq. [10] and then add to it 
the contribution from the cells in the partners list according to Eq. [TTJ 

After these three steps, the only contributions to the magnetostatic field 
that are not counted come from dipole pairs among the leaf cells which are 
near to each other or within the same cell. So for each field point in the leaf 
cell, calculate explicitly Eq. [TBI which sum over all the dipoles in those cells 
in the nearPartners list except the very dipole located at the field point. 

In the majority of implementations of the FMM, the potential is ex- 
panded in terms of spherical harmonics, which complicates the code and 
requires long computing time. Previous attempts were made to expand in 
Cartesian coordinates 14, 15, 0]. Initially, the expansion coefficients were 



calculated either by hand or in a computationally inefficient way, but they 
can in fact be obtained very fast nowadays 0] using Eqs. [12] and [13] The 
disadvantage of Cartesian expansion is that the number of terms below each 
order is larger than that with spherical harmonics, but this inefficiency is less 
than a factor of two up to order 8. We will show later that in most cases in 
the micromagnetic simulation, an order of no more than 6 is sufficient. 

3. Formalism 

In this section, we provide all the necessary formulae to implement the 
above algorithm for nanomagnetic simulation, following Ref . 0] • To simplify 
the formalism, the following shorthand notations are defined: 

n = (n x ,n y ,n z ) 
n = n x + n y + n z 
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n! = njn y \n z \ (1) 
Define the multipole moments of a given charge distribution qf 

^ = ^E^- (2) 

i 

The evaluation of Q n is straightforward for monopoles. For dipoles, how- 
ever, this is not the case. One way to proceed is divide the system always into 
cells which contain only one dipole at the center, based on the fact that only 
Q(i,o,o) ; Q (0,1,0) an d Q (o,o,i) survive for a single dipole m located at the origin, 
and they equal m x , m y and m z respectively. Nevertheless, this method is not 
convenient and even inapplicable for certain systems like ferrofluids. More- 
over, having multiple dipoles in each leaf cell has performance advantages. 
This issue will be discussed in detail in Sec IV. Thus it is important to have 
an easy formula to calculate multipole moments for dipolar system. Previous 
studies have shown complicated formulas for spherical harmonics expansion 
For the multipole moments in Cartesian coordinate, it becomes much 
simpler, i.e. 

Q„ = ^ J^-Vrf. (3) 

i 

In order to derive this formula, we follow the standard procedure to treat 
a dipole rh(f) as a limit of two point monopoles ±q located at r ± d when 
d — ► 0, keeping m = 2qd. After Taylor expansion, Eq. [3] appears. 

If l/\f — fl\ is expanded in Taylor series at rj, the potential V(r) = 
~ ^1 can be expressed in a compact form in terms of the multipole 



moments, 



\f — Tj\ L — ' n! 



£-W> n , (4) 



V(f) = ^D n (r)Q n . (6) 
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In the upward pass (step 2), one needs to shift the origin of the multipole 
moments. The moments Q n about the origin O are related to Q' n about 
position c (with respect to O) by 

= E £tf»-p. ( ? ) 
P 

The Taylor expansion of an arbitrary potential function V(r) is: 

V W = E^ r " ( 8 ) 

n 

Vn = ^V(r)\, =0 . (9) 

In the second downward pass (step 3), one needs to shift the origin of the 
V n when a child inherits V n from its parent. The child's about position c 
(with respect to its parent's origin O) is related to its parent's V n by 

p 

Besides, one also needs to calculate the Taylor expansion coefficient in a 
given cell from its partner's multipole moments. This is achieved by 

K = (-l)"^„ +P (c1Q P , (11) 
p 

where c is the position vector of the center of this cell with respect to the 
center of its partner whose multipole expansion is Q p . It can be shown [?J 
that D n (f) has the form: 

D n(rl = -^T E ^(P) rP ( 12 ) 

p(p=n) 

Thus the problem of calculating D n breaks down to creating the table 
F n (p). Note that by definition of D n , we have: can be calculated by D n . 



D^(f) = -±D. 
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P 

-p x r p -* +2i + (2n + 1 - p x .)r p+x ) (13) 
where x = (1, 0, 0). Meanwhile according to Eq. [12], 

A>+*(0 = ^ E WpV (14) 

p'(p'=n+£) 

Similarly one can obtain the such formulas for -D n +y an d -Dn+z (a per- 
mutation through x, y, z in above equations). Therefore each n, F n (p) can 
clearly be built from F n _p(p) recursively starting from F( 0j0 ,o)((0, 0, 0)) = 1, 
where v is chosen among x, y, z such that n v > 0. Note that for each n, 
F n (p) is nonzero only for a few p's. Thus it is better to store a list of ps 
with nonzero F n (p) for each n and a list of F n (p) value accordingly. Once 
these are set up, D n can be calculated very fast. 

With Eqs. El HI and Q21 we can now discuss the criterion to determine 
near and far cell pairs which is crucial in the first downward pass (step 1). 
Suppose we truncate the terms after the order. The error is of the order 
r Pm /(d — r') Pm 0], where r, r' are the radius of the source and the field cell 
and d is the distance between their centers. In our implementation r = r', 
thus r Pm / (d — r') Pm = [a/ (2 — a)] Pm . Obviously, the series will converge when 
a < 1. Thus a maximum a m < 1 is chosen to determine whether two cells 
are near to each other. In Sec. Ill, we will discuss how to choose it optimally. 

Finally the smooth part of magnetostatic field is given by: 

H s = -W = V ^KVr 11 (15) 
^— ' n! 

n 

And sum over pairs in the nearPartner list C 

H l = H»+ £ ^-3(m,-%)f tJ (ig) 

The total magnetostatic energy E = — m, • Hi in unit of /io/47r. 

4. Performance Optimization 

A benchmark where randomly oriented dipoles are arranged on a simple 
cubic lattice was adopted to test the performance of our program. All the 
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results below are averaged over 50 random configurations and performed 
on an dual Intel P4 3.0GHz with 2G RAM. Fig. ffl demonstrated that the 
complexity of our program is 0(n), and it begins to outperform the brute- 
force calculation around N c = 1000. The small deviations from linearity are 
caused by the change of hierarchy as the size of the systems increases. The 
two parallel lines with different colors and symbols correspond to two sets of 
a m and p m achieving the same accuracy. The shift between them indicates 
that there exists an optimal choice of a m and p m set for a given accuracy. 
In the following analysis, we choose a 32 * 32 * 32 cubic system, but the 
conclusions are independent of system size unless it is too small. 



Figure 1: (color online) Computational complexity for the benchmark system. The black 
square represents the brute- force calculation with complexity 0(n 2 ), whereas the red dot 
and green triangles stand for the FMM with two different sets of parameters (a m , p m ). 
The red line is a least square fit to the data. 

Before discussing of the parameter set, there is another important issue 
which affects the performance significantly and needs to be elucidated first. 
It is the smallest cell size or the average number of particles (s) in the smallest 
cell. The value of s determines the number of levels in the binary tree and 
how many brute-force calculations are required on the finest level. Choosing 
an inappropriate s can deteriorate the performance by a factor of 10. Thus 
it is indeed a nontrivial issue. Obviously, the larger s is, the fewer levels and 
the more brute-force calculations will be required. Because of the fairly large 
overhead of FMM with respect to the brute-force calculation, it is naturally 
expected that s — 1 will not be a good solution. And it is also not good 
to make s too large, for we will lose the power of the FMM then. So there 
must be an optimal number s Q which gives best performance for a given 
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parameter set. Meanwhile, since the overhead increases as the expansion 
order p m increases, the optimal s Q is expected to be a function of p m and 
should increase with p m as well. On the other hand, a m determines the 
number of cells in the "near" region of each specific cell, but the ratio of the 
number of cells in the partners list and the number of cells in the nearPartners 
list should be independent of it. If this ratio remains constant, there is no 
reason to change s Q for fixed p m . Therefore, s Q only depends on p m . Though 
it is not too difficult to come to this conclusion, it is still not clear at all how 
to choose s a . Here, we offer a solution. Table I shows results for a 32 3 system 
with p m = 3, a m = 0.54. In this case, s Q = 8. However as for arbitrary 
system sizes, we cannot divide them freely, the dependence of s on p m is a 
discreet function, and different for different system sizes. We overcome this 
problem by averaging s Q over 30 system size from 30 3 to 59 3 . The result is 
shown in Fig. [2], where each number beside the square points is s Q for p m 
from 2 to 12. The line is a linear fit for those points with p m > 3, indicating 
an approximate power law relationship between s and p m and the power 
is approximately 2. As s Q for a given system size will never be the value 
indicated in the graph anyhow, we reach an empirical equation as follows: 

s = | 0.74/4 p rn > 3 
8 otherwise 

It is very easy to use this equation to estimate s Q . Lastly, for the same 
s Q , one can have different divisions of the system, but it is always better to 
have the smallest cell as close to a cube as possible. 



Tabic 1: Performance for [p m = 3, a m = 0.54) 



s 


time 


1 


34.82 


2 


14.69 


4 


8.28 


8 


6.58 


16 


9.13 


32 


13.69 


64 


20.58 



Based on above discussion, we propose a rule of thumb on how to choose 
the smallest cells: (i) calculate s Q according to Eq. [T71 (ii) the smallest 
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Figure 2: (color online) The optimal number of particles in the smallest cell s as a function 
of expansion order p m . The red line is fitted by the data (p m > 3). The number next to 
each data point is the actual s . 

cells do not have to be cubic and should contain as close to s Q particles as 
possible; (hi) under the above two conditions, the smallest cells should be as 
close to cubic as possible. For example: (i) for 50 3 system with p m = 3, the 
best smallest cell choice is 1.5625*1.5625*3.125 where each cell on average 
contains 7.63 particles; (ii) for 40*40*10 with p m = 4, the best choice is 



2.5*2.5*2.5 (-15.6); (hi) for 32 3 with p m = 5, the best choice is 2*2*4(~16). 



Now that we know how to divide the system, we can optimize the param- 
eter set (p m , a m ). Use e = [a m /(2 — a m )] Pm as the error estimate. Obviously, 
the error decreases as p m increases and a m decreases. To achieve certain 
accuracy, one can have various sets of (p m , a m ) combinations. The computa- 
tional effort, on the other hand, increases as p m increases and a m decreases. 
This claim is easy to understand for p m , since larger p m results in more terms 



Table 2: Performance for e = 0. 5 



Pm «m time 



2 0.36 14.15 

3 0.54 6.71 

4 0.64 6.80 

5 0.71 6.22 

6 0.76 7.36 

7 0.79 8.14 

8 0.81 10.65 
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Figure 3: Optimal computing time versus accuracy. The curve is a guide for the eye. 
The inset table shows the optimal choice of parameter set(p m , a m ) and the corresponding 
computing time for each error e. 



in the expansion. a m , on the other hand, determines how many cells there 
will be in the partners list and nearPartner list. The smaller it is, the more 
cells there are in the lists, resulting in more computing effort. Therefore 
the computing effort varies for different sets of (p m , a m ). Table II shows 
one example for e = 0.05. Obviously, the choice of (p m = 5, a m = 0.71) is 
optimal in this case. The inset table in Fig. [3] gives the optimal choices of 
(Pm, (*m) f° r different accuracies. The curve shows the optimal computing 
time as a function of accuracy. The computing time grows exponentially 
when e < 0.05, so it suggests to choose e = 0.05 if the error is acceptable. 
In fact, the error of magnetostatic energy for e = 0.05 is less than 1% ex- 
cept for configurations whose energy are close to zero. This already reaches 
the requirement of micromagnetics where other uncertainties may make it 
pointless to go beyond such accuracy. Thus (p m = 5, a m = 0.71) is sug- 
gested, if no other requirements exist. Even if higher accuracies are required, 
(.Pm = 6, a m = 0.48) will definitely be sufficient, since the average percentage 
error in this case will be 0.02%. Thus, the inefficiency caused by Cartesian 
coordinates is always less than 2. 

Finally we compare the performance of the Cartesian coordinate FMM 
(CCFMM) with that of the spherical harmonics FMM (SHFMM) and the 
FFT method. For SHFMM, we use the same tree hierarchy as our CCFMM 
so that all the performance difference attributes to the difference of base 
function in expansion: CCFMM uses simple polynomial expansion, while 
SHFMM uses spherical harmonics, which is much more computationally ex- 
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Figure 4: (color online) Comparison of performance of Cartesian coordinate FMM 
(CCFMM) (green solid square: p m = 5, a m = 0.71; red open square: p m = 3, a m = 0.63), 
spherical harmonics FMM (SHFMM) (black dot: p m = 5, a m = 0.71) and FFT method 
(blue star). The inset is a log scale plot. Using Cartesian coordinates makes the FMM 
roughly three times faster. With p rn = 5,a m = 0.71, the performance of CCFMM is 
comparable to the FFT method. When the system is sufficiently large or lower accuracy 
is satisfactory, CCFMM can be superior to FFT. 



pensive. Regarding the error boundary, these two methods is very simi- 
lar as long as we keep the parameter set {s,p m ,a m } the same. They will 
only differ by a small factor of order 1. As to FFT, the routine given in 



Numerical Recipes [16| is used. The results are summarized in Fig. |H As 
expected, the computing time increases linearly for both FMM, and the FFT 
method gives a typical 0(N * log(N)) scaling behavior, which can be seen 
more clearly on the log scale plot. Focusing on the two FMM, CCFMM is 
clearly superior to SHFMM. This is all due to the expensive calculation of 
trigonometric functions of spherical harmonics. Further, cumbersome com- 
plex number calculations are avoided by working in Cartesian coordinates. 
The gain is approximately a factor of three given the chosen parameters. 
This result is comparable to what has been reported in previous literature 
However, differences arise because we use a higher expansion order (p m = 5) 
while earlier authors used (p m = 4). As a rule of thumb: the higher the 
order, the less favorable the cartesian expansion. Another factor that results 
in different speed-ups is the method used to calculate spherical harmonics. 
In our case, we hard-code the expression of spherical harmonics rather than 
using recursion relations. Meanwhile we have also minimized the trigonomet- 
ric functions evaluations. With these optimization of SHFMM though, we 
found the CCFMM clearly exhibits superior performance at relatively lower 
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expansion (p m < 8). To achieve a decent 1% average error, parameter set 
{Pm. = 5, a m = 0.71} is recommended. In this case, the performances of FFT 
and CCFMM are similar. Both the benefit and limitation of FFT is clear. 
It is an exact method, and it is easy to implement. However, FFT requires 
a uniform simple cubit grid and a large padding area for exotic structures 
with open boundary condition. Even though special FFTs can be adapted to 
nonuniform systems, they become problem specific and complicated. FMM, 
on the other hand, is quite generic and natural for nonuniform systems and 
exotic structures with open boundary condition. What's more, it can beat 
FFT when only lower accuracy is needed, as shown by the red open square 
in FigJH 

5. Conclusions 

In conclusion, we have implemented and analyzed the performance of the 
Cartesian coordinate based FMM in dipolar systems. A brief description 
of the algorithm and supplementary equations are provided. In particular, 
we present a simple formula (Eq. [3]) to calculate the multipole moments 
for point dipoles. While carefully analyzing the performance by comparing 
with spherical harmonics FMM and FFT, we have shown that the Cartesian 
coordinate based FMM is appropriate for magnetic simulations due to their 
moderate accuracy requirements. A rule of thumb to decide the optimal 
number of dipoles in the smallest cell is proposed (Eq. IT7I) . and the optimal 
choices of the expansion order and opening angle under different precision 
requirements are discussed. (p m = 5, a m = 0.71) is suggested for nanomag- 
netic simulation generally. As the FMM is a parallel scalable algorithm, it 
can be efficiently implemented in parallel architectures. 

Acknowledgement: We would like to thank Aiichiro Nakano and Yaqi 
Tao for useful discussions. The computing facility was generously provided 
by the University of Southern California high-performance supercomputing 
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A. APPENDIX 

Here we provide the pseudocode for the three recursive function men- 
tioned in Sec. II. Within the framework of object-oriented programming, all 
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these function are public members of an object called "cell". 

void cell::createPartners() 

FOR each cell "A" in the "Partners" list 
IF the cell A is near to this cell THEN 
IF this cell has child THEN 

put the children of A into the "Partners" list of the children of this 
cell; 

erase A from the Partners list of this cell; 
ELSE 

Put A into the NearPartners list of this cell; 
erase A from the Partners list of this cell; 
END IF 
END IF 
END FOR 

IF this cell has child THEN 

put one of its children into the "Partners" list of the other children (do 

this for both children); 

call its children's CreatePartnersQ; 
ELSE 

Add this cell into the NearPartners list of itself; 
END IF 
END FUNCTION 

void cell::updateMoment() 

clear the multipole moments of this cell; 
IF this cell is not the smallest cell Then 
call its children's updateMomentQ; 

sum over its children's multipole moment with shift of origin; 
ELSE IF there is dipoles in this smallest cell THEN 

calculate the multipole moments of this cell directly from the dipole 

distributions; 
END IF 
END FUNCTION 

double cell::updateField() 
IF this cell has parent THEN 

inherit the Taylor expansion from its parent with shift of origin; 
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END IF 

FOR each cell in the "Partners" list 

add to the Taylor expansion of this cell the field generated by A; 
END FOR 

IF this cell has child THEN 
call its children's updateFieldQ; 

calculate the total energy of this cell by the sum of the energy of its 
children; 
ELSE 

calculate smooth local field from the Taylor expansion contributed from 
all far cells; 

FOR each cell "A" in the "nearPartners" list 

add to the local field from each dipoles in A except when the dipole 

is located at the field point; 
END FOR 

calculate the energy of this cell by summing over all the dipole in it and 
divide the energy by 2 to eliminate double inclusion of the interaction 
energy 
END IF 

return the energy of this cell; 
END FUNCTION 
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